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^BSTRACT: One essential question in material sciences is how to bridge the gap between the microscopic 
t^icture and a macroscopic description. The former involves contact forces and deformations, whereas the latter 
-w;oncerns tensorial quantities like the stress or the velocity gradient. 

two-dimensional shear-cell filled with disks of different sizes is examined by means of a "microscopic" 
^iscrete element method (DEM). Applying a consistent averaging formalism, one can obtain scalar- and vector- 
cdields as well as classical tensorial macroscopic quantities like fabric, stress or velocity gradient and, in addition, 
^rjicropolar quantities like curvature or couple-stress. 
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Q INTRODUCTION 

' "the macroscopic balance equations for mass, mo- 
^jpientum and energy can be used for the modeling of 
>the behavior of granular media. However, they rely on 
^Constitutive relations between the physical quantities, 
(^il order to close the system of equations. The deter- 
COnination of both the physical quantities, important for 
•^he behavior of the system, and their inter-relations 
i^io subject of current research (Behringer and Jenk- 
■;;2ns 1997; perrmann et al. 1998| ; [Vermeer et al. 2000| ). 
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)ne possible way to obtain an observable like, for 
t$xample, the stress is to perform experiments with 
,Jf)hoto-elastic m aterial (|Veje et al. 1999 ; Howell and 
^ehringer 1999; Powell et al. 1999D which allows 
QPor non-invasive force measurements. The alternative 
Hs to perform d iscrete particle simulations (C undall 
. ^nd S track 1979; [Hernnann et al. 1998|; L atzel et al. 
^>^000b) and to average over the "microscopic" forces 
^nd contact vectors in the simulation. A typcal snap- 
shot from a two-dimensional DEM simulation is dis- 
played in Fig. |T} 

Besides the trivial definitions for averages over 
scalar and vectorial quantities like density, velocity, 
and particle-spin, one can find differing definitions for 
stress and strain averaging procedures in the literature 
dBathurst and Rothenburg 1988|; Kruyt and Rothen- 
burg 1996; |Liao et al. 1997|; [Calvetti et al. 1997| ; 
f:.atzel et al. 2000bt Pedecker et al. 2OOO0 . 

The aim of this paper is to review recent results 
for tensorial, averaged continuum quantities involving 
also micro-polar aspects ( |Latzel et al. 2UUUa| ; Zervos 
et al. 2000). 



Figure 1: Representative snapshot of the model sys- 
tem. Light and dark grey indicate large and small val- 
ues of the contact potential energy, respectively. 

In section ^ our averaging method is introduced 
and applied to obtain some macroscopic fields in sec- 
tion ^ Section ^ contains the definitions and averag- 
ing strategies for fabric, stress, and elastic deforma- 
tion gradient. The particle rotations are taken into ac- 
count in section |5|, averages are presented for selected 
macroscopic quantities in section ^ and material pa- 
rameters are discussed in section^. 

2 FROM MICRO- TO MACRO-DESCRIPTION 

In a microscopic, discrete picture, the knowledge of 
the forces acting on each particle is sufficient to model 
the dynamics and the statics of the system. Tensorial 
quantities like the stress cr or the strain e are not re- 



quired for a discrete model. In order to establish a cor- 
respondence to continuum theories one has to com- 
pute tensorial fields as well as scalar material prop- 
erties like, e. g., the bulk and shear moduli (C ambou 
et al. 1995; |Kruyt and Rothenburg 1996;_Emeriaul t 
and Chang 1997; |Liao et al. 1997|; |Latzel et al. 2000b|; 



>:.uding et al. 2000ULatzel et al. 2000aD . In the foUow- 
ing, a consistent averaging strategy is introduced. 

2.1 Averaging Strategy 

Most of the measurable quantities in granular materi- 
als vary strongly both in time and on short distances. 
At a fixed point in space, one has either a particle 
or one finds an empty space. Therefore, density flips 
between the two values zero and g^, where the lat- 
ter is the material density. During the "homogeniza- 
tion" presented in the next subsection, we average 
over space and time, in order to reduce the fluctua- 
tions. The averaging over space is an option in sys- 
tems with some translational or rotational invariance 
(simple shear or ring shear apparatus); averaging over 
long times is feasible only in a steady state situation. 
In a ring shear apparatus, as introduced later in sec- 
tion 1^, all points at a certain distance from the origin 
are equivalent to each other and, after several shear 
cycles, the system reaches a quasi-steady state. Aver- 
aging over many snapshots is somehow equivalent to 
an ensemble average. However, we remark that differ- 
ent snapshots are not necessarily independent of each 
other dLatzel et al. 2000b| ) and the duration of a sim- 
ulation might be too short to explore a representative 
part of the phase space. 

The final assumption for the averaging procedure 
presented below, is that all quantities are smeared out 
over one particle. This is a good approximation for 
the density, but not for the stress. However, since we 
aim at macroscopic, averaged quantities, it is not our 
goal to solve for the stress field inside one particle. 
Details of the position dependency inside the particles 
are presumably smeared out due to the averaging in 
space and time. 

2.2 Averaging Formalism 

The mean value of some quantity Q is defined as 

1 



Q = T7 E <y'Q' 



(1) 



with the averaging volume V and the particle volume 
V^. = Q'^ is a quantity attributed to particle 
p, where the quantity Q'^ is attributed to contact c of 
particle p with contacts. The weight Wy accounts 
for the particle's contribution to the average, and cor- 
responds to the fraction of the particle volume that is 
covered by the averaging volume. Since an exact cal- 
culation of the area of a circular particle that lies in 



an arbitrary volume is rather complicated, we assume 
that the boundaries of V are locally straight, i. e. we 
cut the particle in slices, see (|Latzel et al. 2000b|) for 
details. In the following section possible choices for 
observables are summarized and discussed. 

3 MACROSCOPIC FIELDS 

In the following, the averaging formalism in Eq. (|T]) 
is applied to obtain various macroscopic quantities. 
For example, the quantity Q'^ can be chosen as in 
order to obtain the density, as (fv^ in order to obtain 
the momentum density, or as {l/2)(f(v'''y for the 
kinetic energy density. In table [I], the macroscopic 
fields as well as particle-attributed quantities are 
collected. In this study, two tensors with no symbol in 
between mean a dyadic product, whereas '-'is used 
for the scalar product, i. e. the order-reduction by one 
for each of the two tensors at left and right. 

Table 1 : Macroscopic fields computed using the aver- 
aging formalism in Eq. (P using particle properties 



macroscopic quantity 




volume fraction u 


1 


density g 




flux density uv 




momentum density gv 


gPvP 


kinetic energy density Igv'^ 
dynamic stress tensor cr 


yp{vpy 


gPv^vP 



3 . 1 Volume Fraction 

As a first example for an averaged scalar quantity, the 
local volume fraction u, see Tab. [T], is directly related 
to the local density via g ^ (fy . 

3.2 Flux Density 

As a tensorial quantity of first order (vector), the ve- 
locity field V is obtained when dividing the flux den- 
sity uv by the volume fraction u. The momentum den- 
sity is thus gv gPuv . 

3.3 Velocity Gradient 

Given the velocity field, it is possible to compute the 
velocity gradient Vv, a tensor of second order, by 
means of numerical differentiation. With Vi> one has 
the components of the symmetric deformation rate 
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(2) 



and the anti- symmetric continuum rotation rate 



W. 



af3 



d(3 da 



(3) 



with the coordinates a and In the special case of 
cylindrical symmetry with radial direction r and an- 
gular direction 0, Eqs. (§) and (||) reduce to 



and 
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r 



dr r 



(4) 



(5) 



see ( [Luding et al. 2OO0|) for details and (Z ervos et al. 
2000) for a similar approach. 

3.4 Energy Density 

The kinetic energy density is only a function of the 
velocity of the particles in the averaging volume and 
it is related to the trace of the dynamic stress tensor. 
However, we do not discuss this quantity here, since 
the system examined later is quasi-static so that the 
dynamic cont ributions are of minor importance (L ud- 
ing et al. 2000; |Latzel et al. 2000a^ 



4 MACROSCOPIC TENSORIAL QUANTITIES 
In this section, the averaged, macroscopic tensorial 
quantities in our model system are introduced. The 
fabric tensor describes the statistics of the contact di- 
rections, the stress tensor describes the stress distri- 
bution due to the contact forces, and the elastic defor- 
mation gradient is a measure for the elastic, reversible 
deformations due to the stress. 

Table 2: Macroscopic tensorial quantities computed 
with Eq. (jl]) using contact-attributed properties pre- 
averaged over single particles 



macroscopic quantity 


QP 


fabric tensor F 
static stress tensor a 
deformation gradient e 


c=l 
1 CP 

pcipc 

c=l 

h '-'^ 

— V A^n'p^ ■ F-^ 



4. 1 Fabric Tensor 

In assemblies of grains, the forces are transmitted 
from one particle to the next only at the contacts of 
the particles. Therefore, the local geometry and direc- 
tion of each contact is important ( p^atzel et al. 2000b[ ). 
The fabric tensor, in our definition, involves the con- 
tact normal vector n'^, related to the so-called branch 
vector via = apin^, with particle radius a^. For each 
contact the dyadic product n^n^ is used, i.e. a de- 
generate tensor of order two with non-zero value in 



n'^-direction only. In average over many contacts (and 
over many particles) one obtains the fabric tensor. 

The fabric tensor in table |^ is symmetric by def- 
inition and thus consists of up to three independent 
scalar quantities in two dimensions. The first of them, 
the trace (or volumetric part) Fy = tr(F) = (F^ax + 
F-ami), is the contact number density, with the major 
and the minor eigenvalues F^ax and F^in, respectively. 
With other words, one obtains the relation tr(F) = uC 
with reasonable accuracy for monodisperse particles, 
where C is the average number of contacts per par- 
ticle. For polydisperse size distribution functions, a 
multiplicative correction factor was recently proposed 
dMadadi et al. 20011 ), which depends on the first three 
moments of the size distribution function only. 

The second scalar, the deviator Fo = F^ax — F^m, 
accounts for the anisotropy of the contact network to 
first order, and the third, the angle (pp, gives the ori- 
entation of the "major eigenvector", i. e. the eigenvec- 
tor corresponding to F^ax, defined with respect to the 
radial outwards direction for cylindrical symmetry. 
In other words, the contact probability distribution is 
proportional to the function Fy + F^ cos(2(</) — (pp)) 
dLatzel et al. 2000b| ), when averaged over many par- 
ticles. Note that this approximation, a fourier expan- 
sion up to second order, is not always sufficient, i. e. a 
more general fabric tensor of rank four should be used 
dMehrabadi et al. 1988| ; perrmann et al. 1998| ). 



4.2 Stress Tensor 

The stress tensor cr is a measure for the force f = cr ■ 
n, acting on an imaginary surface with normal n, per 
unit area. A force is equivalent to some momentum 
transfer per unit time, so that also the momentum flux 
density contributes to the stress in a dynamic way. The 
dynamic stress, see table [1], corresponds to the stress 
in an ideal gas, whereas the static stress, see table 0, is 
caused by the transfer of forces from a contact to the 
center of mass of a particle. 

In order to account for both the stress and the dis- 
tance of transfer, the static component of the stress 
tensor ( piothenburg and Selvadurai 1981 ; Satake and 



Jenkins 1988; Kxuyt and Rothenburg 1996|; L atzel 
et al. 2000b) is defined as the dyadic product of the 
force acting at contact c with the corresponding 
branch vector, see table 0. Note that the formerly in- 
troduced surface quantity is here expressed as the sum 
over all particles in the averaging volume, with the re- 
spective weight factor in Eq. ([H), for a detailed deriva- 
tion see dLatzel et al. 2000b| ). 



4.3 Elastic Deformation Gradient 

In order to reach the final goal, i. e. to measure the 
material properties of a granular ensemble, one is in- 
terested, e. g., in the stress-strain relationship of the 
material. The strain e can be obtained by time Integra- 



tion of the velocity gradient, see subsection and 
subsequent symmetrization and linearization. Here, 
an alternative quantity is introduced by application 
of "Voigt's hypothesis", i. e. assuming that the defor- 
mation is uniform and that every particle displace- 
ment conforms to the corresponding mean displace- 
ment field, but fluctuates about it ( |Liao et al. 1997| ; 
Latzel et al. 2UUUbp . The deformation gradient in table 
g|is calculated for the special case of two-dimensional 
disks with height h. 

This relates the actual deformations to a virtual, 
stress-free reference state where all contacts start to 
form, i. e. t he particles are just touching, see (L atzel 
et al. 2000b) for a detailed derivation. The result is a 
non-symmetric tensor e, which is not the strain. In- 
stead, we refer to it as the elastic deformation gradi- 
ent, since it accounts only for reversible (elastic) de- 
formations. 



5 ROTATIONAL DEGREES OF FREEDOM 
Due to the particles' surface roughness, forces are 
transmitted also tangentially and thus granular parti- 
cles will rotate. Therefore, also micro-polar macro- 
scopic quantities related to the rotational degrees of 
freedom a re of interest ([Latzel et al. 2000a|; Z ervos 
et al. 2000). The particle rotation, the couple-stress 
and the curvature are discussed in this section. 

Table 3: Tensorial quantities, connected to the rota- 
tional degrees of freedom, computed with Eq. ((T|) 



macroscopic quantity 




spin density z/cj 
couple stress M 

elastic curvature k 


1 CP 
c=l 

7lh 

c=l 



5 . 1 Particle Rotation 

The mean angular velocity of the particles is obtained 
by using Qp = ujP, see table ||. Note that w also con- 
tains the continuum angular velocity 



a;, = I X W 



(6) 



i.e. the vector product ' x ' of the continuum rotation 
rate from Eq. (^ and the unit tensor I. The excess 
rotation, or particle eigen rotation, with respect to the 
underlying mean motion, is thus 



(jj = (jj — (jj 



(7) 



5.2 Couple stress 

In the fra mework of a Cosserat continuum (Z ervos 
et al. 2000) one has, in addition to the stress, also the 
couple stress M. When an applied stress leads to a de- 
formation, an applied couple-stress causes a rotational 
motion. It can thus be defined in analogy to the stress 
by replacing the force by the torque due to the tan- 
gential component of the force, see table ^ In a two 
dimensional system with cylindrical symmetry, only 
the components Af^^ and M^,^ of the tensor survive. 
Note that M = 0, when the sum of the torques acting 
on spherical particles vanishes by definition in static 
equilibrium. 

5.3 Elastic Curvature 

In analogy to the elastic deformation gradient, we de- 
fine the curvature k by replacing the overlap with the 
respective axial vector, x A^'^, see table ^ This 
leads to a measure for the reversible or "frozen in" 
rotations in the system. Consequently, in static equi- 
librium of spherical particles, one has k = 0. 

6 RESULTS 

The elementary units of granular materials are meso- 
scopic grai ns. With DEM ([Latzel et al. 2DDDE|; L atzel 
et al. 2000a) the grains are treated as rigid particles 
but their local deformation at their contact points is 
realized as virtual overlaps. We relate the interaction 
forces to this overlap 5 and to the tangential displace- 
ment of two particles during contact. The force laws 
used are material dependent and have to be validated 
by comparison w ith experimental measurements (F o- 

Falcon et al. 



i. e. the total mean spin minus the continuum spin. 



erste r et al. 1994; [La bous et al. 1997 
1998; [Luding 1998[ ). 



6.1 The Model System 

In the simulations presented in this study, a two- 
dimensional Couette shear-cell is used. N disks are 
confined between an outer ring and an inner ring 
with radius Ro and respectively. The particles 
are of slightly different size rfsmaii = 7.42 mm and 
enlarge = 8.99 mm, in order to reduce ordering effects. 
These boundary conditions are based on an experi- 
ment ( [Howell and Behringer 1999[ ; IVeje et al. 1999| ); 
for more d etails and other simulations, see (L atzel 
et al. 2000b; [Luding et al. 2000^ 

The outer wall is fixed and the inner wall rotates 
and thus introduces a slow shear deformation to the 
system. The simulations are started in a dilute state 
with an extended outer ring while the inner ring al- 
ready rotates counter-clockwise with constant angular 
velocity Vt = 2n/Ti = 0.1 s^^ and period = 62.83 s. 
The radius of the outer ring is reduced within about 
two seconds to reach its desired value Ro and there- 
after it is kept fixed. Averages are performed after 
about three rotations at t = 180 s (to get rid of the 



arbitrary initial configuration), and during about one 
rotation, until t = 239 s. In Fig. [I], a typical snapshot is 
displayed, where the potential energy of the particles 
is coded in greyscale so that the stress-chains become 
visible. 

Different global volume fractions z/ were examined 
in the simulations. Here, we present data from three 
different simulations A, B, and C with u = 0.8084, 
0.8149 and 0.8194, respectively. For the different sim- 
ulations the number of large and small particles is 
N,^^n- 2511, 2545, 2555 and iVi^.ge: 400, 394, 399. 
For the calculation of the global volume fraction, the 
small particles glued to the wall are counted with half 
their volume only, and thus contribute with z/waii = 
0.0047 to u. For more details and material parame- 
ters, see dLatzel et al. 2000bD . 



6.2 Volume Fraction 

In Fig. 1^, the volume fraction is plotted against the 
rescaled distance from the inner ring. Starting from 
an initially uniform volume fraction, a dilated shear- 
zone forms near to the inner wheel as a consequence 
of the applied shear. This effect is less pronounced 
for higher initial global densities and, in the outer re- 
gion of the shear-cell (f > 10), the structure of the 
packing remains frozen, i. e. not much reorganization 
takes place within the duration of the simulation. 
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Figure 2: Volume fraction v, plotted against the di- 
mensionless distance from the origin f = (r — Ri) / d, 
with the typical diameter d = 8 mm, for different ini- 
tial global densities z/ 



6.3 Velocity Field 

In Fig. |3|, the tangential velocity is plotted against the 
rescaled distance from the inner ring. The simulation 
data are fitted by a function f</,(r) = foexp(— f/s) 
with Vo-. 0.670, 0.756, 0.788 and s: 1.662, 1.584, 
1.191, thus showing an exponential profile corre- 
sponding to the shear band. The shear band, has a 
width of a few particle diameters, before the velocity 
f reaches the noise level. 
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Figure 3: Tangential velocity normalized by the 
velocity of the inner wheel ^IRi, plotted against f. 
The lines are the fits to the simulation data for 
f = 0.25 to 8.1, using v^{r) = voexp{—f/s) with Vo'. 
0.670, 0.756, 0.788 and s: 1.662, 1.584, 1.191, for in- 
creasing density, as given in the inset, respectively 



6.4 Fabric Tensor 

The trace of the fabric tensor is a measure for the con- 
tact number density, whereas the deviator of the fabric 
is a measure for the anisotropy of the contact network 
dLatzel et al. 2UUUb| ; [Latzel et al. 2l}DDaD , see Fig. |. In 
our situation, the contact number density is reduced in 
the shear zone and the anisotropy, i. e. the deviatoric 
fraction, is increased in the shear zone, but remains 
below 20 per-cent. 
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Figure 4: Components of the fabric tensor tr(F) 
(isotropic) and dev(F)/tr(F) (anisotropic) for differ- 
ent numbers of binning intervals 



6.5 Stress Tensor 

In Fig. ^ the static stress components are plotted. 
Here, the diagonal elements are almost constant, 
whereas the off-diagonal elements decay proportional 
to r~^, as indicated by the lines, in consistency with 
continuum theoretical equilibrium conditions (L uding 
et al. 2000). The trace of the stress tenso^.e. the 
volumetric stress, is almost constant over the whole 
shear-cell besides fluctuations. In contrast, the devia- 
toric fraction ao/ay decays with increasing distance 
r from the inner ring, similar to the behavior of the 
deviatoric fraction of F. 
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Figure 5: Components of the static stress <t plotted 
against the distance from the center r in log-log rep- 
resentation. 

6.6 Tensor Orientations 

The orientations of the tensors F, cr, and e (data not 
shown here) are tilted from the radial outwards direc- 
tion into s hear direction (^atzel et al. 20Q0b|; L uding 
et al. 2000; [Latzel et al. 2000aD . Interestingly, the ori- 
entations of the three tensors are different, indicating 
an anisotropic material. All orientation angles show 
the same qualitative behavior, however, the fabric is 
tilted more than the stress which, in turn, is tilted more 
than the deformation gradient. Thus, the three tenso- 
rial quantities examined so far are not co-linear. 

6.7 Angular velocity field 

In Fig. H the (^-components of the) macroscopic par- 
ticle rotations uj, the continuum rotation uj^ and the 
particle excess-rotation uj* = uj — uJc, are displayed. 
Both the particle- and the continuum-angular veloc- 
ity decay exponentially with increasing f, paralleling 
the behavior of the velocity v^. The inset of Fig. § 
shows an oscillation of the excess-rotation near the 
inner wheel, from one disk layer to the next. This is 
due to the fact that the disks in adjacent layers are able 
to roll over each other in the shear zone. 

7 MATERIAL PARAMETERS 
7.1 Bulk Stiffness 

The material stiffness E = tr(cr)/tr(e) is here defined 
as the ratio of the volumetric parts of stress and elas- 
tic deformation gradient. In a mean field estimate it is 
found to be proportional to the "microscopic" con- 
tact stiffness k'^ and the trace of the fabric tensor: 
E oc (A;y27r)tr(F). In Fig. [7] the rescaled stiffness 
of the granular material is plotted against the trace 
of the fabric. Note that all data collapse almost on a 
line, but the mean-field value underestimates the sim- 
ulation data by some per-cent. The few data points 
which deviate most are close to the boundaries. The 
deviation from the mean field prediction (solid line in 
Fig. ^ seems to disappear in the absence of shear. 
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Figure 6: Angular velocities u (solid line) and cUc 
(symbols) of the particles and the continuum, plotted 
against the scaled radial distance (from simulation B). 
The dotted line is cOc as obtained from the fit to v^, 
see subsection In the inset, the excess-spin is dis- 
played for all simulations 
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Figure 7: Scaled granular bulk stiffness 2'nE/k[^, plot- 
ted against tr(F). Every point corresponds to one of 
150 rings dividing the system 

7.2 Shear Stiffness 

The typical shear stiffness of the material is de- 
fined as the ratio of the deviatoric parts of stress 
and elastic deformation gradient G = dev(cr)/dev(e) 
and scales in a crude mean field approximation as 
G oc (A;^/7r)tr(F). In Fig. |]the ratio of the deviatoric 
parts of stress and strain is plotted against the trace of 
the fabric. Like the bulk stiffness, both quantities are 
proportional, at least for points near or within the di- 
lute shear band. In the denser outer part of the shear- 
cell, the particles are strongly inter-locked and thus 
resist much more against shear, so that G diverges. 
The critical contact number density grows with in- 
creasing global density, i. e. with increasing stress. 

7.3 Torque Resistance 

Since we are interested in the role the rotational de- 
gree of freedom plays for the constitutive equations. 




Figure 8: Scaled granular shear resistance -nG/k'^ = 
dev(cr)/dev(e) plotted against tr(F) from all simula- 
tions. The line indicates the identity curve 
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Figure 9: Torque resistance plotted against f 



we define a new material parameter, the "torque resis- 
tance" fic, as the ratio of the magnitudes of the cou- 
ple stress and the curvature components. This quan- 
tity describes how strongly the material resists against 
applied torques. In Fig. ^ is plotted against the 
scaled distance from the inner ring. In the dilute re- 
gions near to the shear zone, where the particles are 
able to rotate more easily, fi^ is smaller than in the 
denser regions, where the particles are interlocked and 
thus frustrated. This is consistent with the results for 
increasing global densities, i. e. the torque resistance 
increases with density. Note that the strongest fluctua- 
tions are due to the division by small values and have 
no physical relevance in our interpretation. 

8 SUMMARY AND CONCLUSION 

A micro-macro averaging procedure was presented 
and specific boundary conditions were chosen to al- 
low for averaging over large volumes and long times. 
A shear band is obtained, localized close to the in- 
ner, rotating cylindrical wall. The local configurations 
changed rather rapidly in the shear band, whereas the 
system is almost frozen in the outer part. 

The averaging strategy used assumes the quantities 



to be homogeneously smeared out over the whole par- 
ticle which is cut in slices by the averaging volumes. 
The material density, i. e. the volume fraction, the ve- 
locity field, the fabric tensor, the stress tensor and the 
elastic, reversible deformation gradient were obtained 
from DEM simulations. The volumetric stress is con- 
stant in radial direction while its deviatoric fraction 
and also the mean velocity gradient decay with in- 
creasing distance from the inner wall. The ratio of the 
volumetric parts of stress and strain gives the material 
stiffness of the granulate, which is small in the shear 
band, due to dilation, and larger outside. 

In the shear-cell, large deviators, i. e. anisotropy, 
of all tensorial quantities is evidenced. All tensors 
are tilted from the radial direction. The system or- 
ganizes itself such that more contacts are created to 
act against the shear and also the shear resistance in- 
creases with the contact density. An essential result 
is that the macroscopic tensors are not co-linear, i. e. 
their orientations are different. Thus, the material can- 
not be described by a simple elastic model involving 
only the two Lame constants (or bulk modulus and 
Poisson's ratio). 

Finally, the particles' angular velocity was aver- 
aged in analogy to the particle velocity. Subtraction 
of the continuum rotation from the particle rotation 
leads to the excess-eigen rotation of the particles with 
respect to the mean rotation, in the spirit of a micro- 
polar or Cosserat continuum theory. In analogy to 
the stress and elastic deformation gradient, we de- 
fined couple stress and curvature. The quotient of 
the respective non-zero components gives a "torque- 
resistance" which increases with increasing local den- 
sity and stress. 

Current research concerns the application of the av- 
eraging formalism to other boundary conditions, the 
use of more realistic particle interaction models, and 
the measurement of other macroscopic quantities not 
discussed here. 
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